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ABSTRACT 

We present the first fully relativistic investigation of matter fallback in a supernova. We 
investigate spherically symmetric supernova fallback using a relativistic radiation hydrodynamics 
' Lagrangian code that handles radiation transport in all regimes. Our goal is to answer the 

fundamental question: did SN1987A produce a black hole and, if so, when will the hole become 
detectable ? We compute the light curve, assuming that a black hole has been formed during the 
explosion, and compare it with the observations. Our preliminary calculations lack radioactive 
energy input and adopt a very simple chemical composition (pure hydrogen). As a result, our 
computed models cannot fit the observed data of SN1987A in detail. Nevertheless, we can show 
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that, during the first hours, the accretion flow is self-regulated and the accretion luminosity 
stays very close to the Eddington limit. The light curve is completely dominated, during the first 
few weeks, by the emission of the stellar envelope thermal energy, and resembles that obtained 
in "standard" supernova theory. Only long after hydrogen recombination takes place is there 
even a chance to actually detect radiation from the accreting black hole above the emission of 
the expanding envelope. The presence of a black hole is thus not inconsistent with observations 
to date. Because of the exponential decay of the 44 Ti radioactive heating rate, the date of the 
\ emergence of the black hole is not very sensitive to the actual parameters of the models and 

turns out to be about 1000 years. The bulk of the emission then is expected to be in the visible 
band, but will be unobservable with present instrumentation. We discuss the implications of our 
\ results in connection with the possible emergence of a black hole in other supernovae. 

CO 

Subject headings: Accretion, accretion disks — hydrodynamics — methods: numerical 
radiative transfer — relativity — supernovae: individual (SN1987A) 

1. Introduction 

On the 23rd of February 1987, Shelton and Jones announced the discovery of a supernova in the 
Large Magellanic Cloud, SN1987A (Kunkel et al. 1987). Historically, this was the brightest supernova 
observed after that recorded by Kepler in 1604 (SN1604). It was also the first supernova to be observed in 
every band of the electromagnetic spectrum and the first detected through its initial burst of neutrinos. 
Its relatively close distance (~ 50 Kpc) has offered a unique opportunity to observe a supernova in great 
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detail and with a variety of detection techniques. [For extensive reviews on SN1987A see e.g. Arnett et al. 
(1989), McCray (1993) and references therein.] Since its appearance, SN1987A has confirmed many aspects 
of the theory of type II supernovae (see e.g. Dar 1997). As predicted, most of the energy (~ 10 53 erg) was 
released in ~ 10 s by the cooling of the protoneutron star and was radiated in form of neutrinos, detected 
by the Mont Blanc, Kamiokandc, 1MB and Baksan underground detectors. The strong shock produced by 
the core bounce heated and pushed outwards the stellar envelope, which subsequently emitted its internal 
energy as photons: the integrated light emission was <~ 10 49 erg, while the kinetic energy of the expanding 
shells was <~ 10 51 erg. The light curve can be fitted quite well by the "standard" type II supernova theory 
(see e.g. Woosley 1988; Arnett 1996). After the emergence of the shock wave, the emission is dominated 
by the diffusion luminosity of the expanding stellar envelope during the first 30-40 days. At <~ 40 days, 
the hydrogen envelope starts to recombine and a large amount of internal energy is released by the inward 
motion of the recombination front. After most of the envelope has recombined at day ^160, the light curve 
starts to be dominated by the radioactive decay of heavy elements synthesized during the explosion. Some 
predictions of the "standard" theory were not confirmed and a number of key issues remain to be answered, 
such as the nature of the supernova progenitor, the explosion mechanism and, perhaps most important, the 
nature of the compact remnant left over after the explosion. 

Stellar evolution calculations show that stars with a main sequence mass in the range 8-19 Mq finish 
their lives with a compact core of <~ 1.4M Q . As shown by Woosley & Weaver (1995), for these stars the 
amount of material that falls back toward the core in the aftermath of its collapse is negligible, so it is likely 
that they give birth to neutron stars. Stars with main sequence mass larger than M c ~ 25-30 Mq have 
more massive cores (~ 2M Q ) and undergo accretion of a significant amount of matter, so that the remnant 
mass left over after the explosion is larger than 3M , probably leading to the formation of a black hole. 
The value of M c is somewhat uncertain and depends quite sensitively on the mass loss prior to the explosion 
(reduced core mass) and on the explosion energy (amount of fall back). The fate of stars in the intermediate 
range 19-25 Mq is far less obvious. At the end of their evolution, they have core masses around 1.6-1.8 
Mq and a variable amount of matter, 0.1-0.3 Mq may fall back because of the hydrodynamic interaction 
of the outgoing supernova shock with the expanding envelope (Woosley & Weaver 1995). Typically, these 
stars leave a central compact remnant of 1.7-2.1 Mq (baryonic mass). Then their fate depends critically on 
the equation of state at nuclear matter densities that fixes the maximum mass M cr i t above which no stable 
neutron star configuration can exist. Adopting the parametrization of the nuclear matter equation of state 
by Lattimer & Swesty (1991), M crit turns out to be M crit <J 2.5M Q (gravitational mass). In this case, it 
would be unlikely that these stars could form a black hole. However, Thorsson, Prakash & Lattimer (1994) 
have proposed a much softer equation of state as a consequence of K~ condensation. They find a maximum 
mass M cr it — 1-5M Q (gravitational mass). Recently, modern models of nuclear forces and accurate Monte 
Carlo modeling of nuclcon interactions indicate that M crit = 1.8-2.2 Mq (Pandharipandc 1997). Based 
on the findings of Thorsson, Prakash & Lattimer (1994), Brown & Bethe (1994) and Woosley & Timmes 
(1996) have proposed different scenarios that may lead to the formation of a black hole after the explosion 
of a 19-25 Mq star. According to Brown & Bethe (1994), after core bounce a protoneutron star would form 
that could remain stable for about ~ 12 s (just the time necessary to release its internal energy in form of 
neutrinos) and then would collapse to a black hole, as early suggested by Wilson et el. (1986) and Woosley 
& Weaver (1986). Woosley & Timmes (1996) point out that the formation of a black hole could be driven 
by the matter which falls back in the first few hours after the explosion, as already emphasized by Colgate 
(1971, 1988) and Chevalier (1989). 

Interestingly, the mass of the progenitor of SN1987A (M = 18-21 Mq) falls in the range of masses 
where the outcome of core collapse is uncertain. So, from a theoretical point of view, we do not know if a 
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neutron star or a black hole has been formed during the explosion. The main uncertainties are the value 
of the critical mass above which no stable neutron star can exist and the amount of matter that falls back 
during the explosion. If there was no fallback in SN1987A, the emission from a radio pulsar most probably 
would not be detectable directly for a very long time, unless there are holes in the outflow. Observations 
also rule out the presence of any optical pulsar. The energy deposition from a pulsar could alter the light 
curve, though; observations to date limit any energy input from a pulsar to be lower than the contribution 
from radioactive decay (~ 10 36 erg s" 1 , implying B/P 2 ~ 10 14 Gs~ 2 for magnetic dipole emission by a 
pulsar with field strength B and rotation period P). However, high energy X-rays from a pulsar radiating 
a hard spectrum should have become observable in a few months (McCray, Shull and Sutherland 1987) 
although the 1-3 keV thermal X-ray emission from the cooling of a neutron star would be attenuated for 
about 100 yr (Chernoff, Shapiro and Wasserman 1989). Such hard X-rays were never observed. 

If a compact object of any kind is present, it is probably accreting from the progenitor stellar material, 
since the inner part of the ejecta have velocity smaller than the escape velocity. During the first phases 
after the explosion, the accretion rate should be as high as 10 8 in units of the Eddington accretion rate 
(Chevalier 1989). On the other hand, as mentioned above, the bolometric light curve observed to date (ten 
years after the explosion) can be explained by the "standard" theory of Type II supernovae and does not 
show any feature related to the possible presence of a central accreting compact object. Using analytic and 
numerical methods, Chevalier (1989) and Houck & Chevalier (1991) have investigated the problem of fall 
back onto a neutron star and found that, for radiation pressure -limited accretion, the luminosity emitted 
after 3-4 years should approach the Eddington value, L^dd ~ 10 38 erg s _1 for gas opacity dominated by 
electron (Thomson) scattering. Since at <~ 1500 days the bolometric luminosity of SN1987A was ~ 10 37 erg 
s _1 , they argue that either the accreting envelope becomes dynamically unstable or SN1987A contains a 
black hole. On the other hand, Chen & Colgate (1995) noted that, at the typical temperature and density 
of the ejecta after few years, the opacity in the accreting gas enriched of heavy elements is much larger 
than the Thomson opacity. So, they conclude that a neutron star accreting at the actual (non-Thomson) 
Eddington limit may be present. Furthermore, if the value of the critical mass falls in the range 1.8-2.2 
Mq, as seems to be indicated by recent developments in the nuclear many-body theory (Pandharipande 
1997), it appears that the amount of fallback is a key element in determining whether or not a stellar black 
hole has formed. SN1987A may or may not harbor a black hole but observations alone have not resolved 
this issue yet. It is therefore important to determine what extra luminosity an accreting central component 
would produce and when its presence might be discernible. The amount of fallback depends on the main 
parameters of the explosion while the visibility of the remnant on the detailed hydrodynamical evolution of 
the flow and of the radiation field. 

This is the first of a series of papers in which we address a number of relevant questions: can we infer 
from the light curve whether a stellar black hole has formed in the aftermath of a supernova explosion ? Did 
SN1987A produce a neutron star or a black hole ? We will consider these issues by assuming that either a 
black hole or a neutron star has formed and then will compute the resulting light curve varying the main 
parameters of the explosion. Our focus will be to isolate the contribution from the central compact object. 

In this paper we present the first self-consistent, fully relativistic investigation of supernova fallback in 
presence of a black hole. It is the relativistic generalization of the work by Colpi, Shapiro & Wasserman 
(1996, hereafter CSW) in which the hydrodynamical properties of a fluid accreting onto a central remnant 
from an initially expanding cloud were first explored. Here, we include self-consistently the transfer of 
radiation. Due to the complexity of the problem we will not attempt to investigate a model with realistic 
chemical composition (that will be presented in a forthcoming paper) but, instead, we will consider a 
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number of simpler preliminary scenarios to guide our understanding. We will show that, during the first 
hours, the accretion flow is self-regulated and the accretion luminosity stays very close to the Eddington 
limit. However, during the first few weeks the light curve is completely dominated by the emission of 
the stellar envelope internal energy, giving a light curve that resembles that obtained in the "standard" 
supernova theory. Although in this paper we apply our calculation specifically to SN1987A, our results have 
general validity and can be used to study the light curve and the observational signatures of black holes in 
other supernovae. 

We have investigated spherically symmetric supernova fallback using a general-relativistic, radiation 
hydrodynamic Lagrangian code. Our code can handle the transfer of radiation from the first phases 
immediately after the supernova explosion (when photons diffuse through a high temperature, expanding 
cloud) to the late evolutionary stages (when the hydrogen envelope has recombined and most of the ejecta 
are transparent). The main challenge we have confronted is the enormous dynamic range in the relevant 
physical timescales entering the problem. To this end, we have implemented a multiple timestep procedure 
that allows us to integrate different radial regions at different rates (a primitive "temporal adaptive mesh"). 
This technique enables us to speed up the calculation by almost an order of magnitude. 

The plan of the paper is the following: in Section |2| we present the general relativistic equations of 
radiation hydrodynamics used in the calculation. The finite difference form of these equations is presented 
in the Appendix. Section ^ describes the numerical method used to solve the equations of radiation 
hydrodynamics and the boundary conditions adopted. In Section |4| we present in detail the multiple 
timestep procedure employed in our calculation. Section ^ illustrates the set up of the initial conditions. 
In Section ^ we present the numerical results for various simplified examples and test cases. Section ^ is 
devoted to the semi-analytic calculation of the late-time light curve. Finally Section || contains a discussion 
of our preliminary results in connection with SN1987A and with the diagnosis of black holes in supernova 
explosions. 



In this section we present the equations of relativistic radiation hydrodynamics in spherical symmetry 
for a self-gravitating matter fluid which is interacting with radiation [for a review of the derivation of 
these equations see Zampieri 1995 and Zampieri, Miller & Turolla 1996]. All the equations will be written 
in the frame comoving with the fluid flow. The 4-velocity of an element of fluid u a will be evaluated in 
the Eulerian frame, defined as the reference frame at rest with respect to the background (spherically 
symmetric) metric. In this frame, the areal radius r is taken to be the independent radial (Schwarzschild) 
coordinate. We now introduce the spherically symmetric, comoving-frame line element 



where t and fi are the Lagrangian time and the comoving radial coordinate (taken to be the rest mass 
contained within a comoving spherical shell) and a and b are two functions of t and fi which need to be 
computed from the Einstein field equations. Here and throughout we adopt geometrized units and set 
c = G = 1. Using the line element ([!]), the complete system of radiation hydrodynamics equations in the 
frame comoving with the flow along with the Einstein field equations can be cast into the form (Rezzolla & 
Miller 1994; Zampieri, Miller & Turolla 1996) 



2. Equations 
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where 



M„ = 47rr 2 r iA1 (e + w + ^u>i) (7) 

U= ~a~ < 8 > 
is the radial component of the fluid 4-velocity measured in the Eulerian frame, T = r^/b = (1+u 2 —2M / r) 1 / 2 , 
M represents the effective gravitational mass-energy (for black hole + gas + radiation) contained within 
radius r and p, e, p and h = (e +p)/p are the rest-mass density, total mass-energy density, pressure 
and enthalpy of the gas flow, respectively, as measured in the comoving frame. Here sq and si are the 
radial PSTF moments of the source function that account for energy and momentum exchange between 
the fluid flow and the radiation field. Their actual forms depend on the radiative processes considered (see 
below). Finally, wq, Wi and u>2 denote the first three radial PSTF moments of the specific intensity / of the 
radiation field, measured in the comoving frame and given by 

kH2k + 1) f 

= 27T .'\ ; ' / IP k {p)dfl, (9) 



(2fc + l)!! 

where p, is the cosine of the angle between the photon propagation and radial directions and P k (fi) is the 
Legendre polynomial of order k. In terms of the "classical" moments of the specific intensity J, H and K, 
we have: Wq = Air J, w± — AirH, u> 2 = Ait(K — J/3). 

The radiation hydrodynamic equations must be supplemented with the radiation moment equations. 
In spherical symmetry with the line-element (|l]), the first two frequency-integrated PSTF moments of the 
relativistic transfer equation can be written (Zampieri, Miller & Turolla 1996) 

1 / .4/3 8/s\ , 1 / 2 2\ , / bt r,t 



64/3 r 8/3 



(u^/V/ 3 ) ( + J-j^oVXp + (y - o-o = , (10) 

^2 M V X t + ^b^ ^^ + ^3 ( W2ar3 )^ - asi = . (11) 

In equations ( fTpf ) and ( pT| ) wq and w\ have the dimensions of energy density. To close the moment 
equations, we employ interpolations of the Eddington factors for a spherically symmetric, static atmosphere 
in radiative equilibrium, that fit the data from Hummer & Rybicki (1971) to within 30%. The expression 
for the Eddington factor / = W2/W0 and the Eddington boundary factor g = wi(r out ) /wo(r ou t) are given by 



l-0.59(^- 

3 \r out 



' (12) 



l + 3r ' 



= 1-0.423 (-2-) , (13) 
\r ut ) 

where r is the optical depth (see below), r± is the radius at which r = 1 and r out is the outer radial 
boundary of the flow. Note that < / < 2/3 — 0.59(ri/r onf ). The use of equations (O) and (G~3h is not 
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exact, especially for a moving gas. However, adopting any kind of iterative procedure to calculate the 
Eddington factors precisely (see e.g. Shapiro 1996) would be too expensive from the computational point 
of view. Moreover, the error introduced by equations ( fl2|) and ( |l3| ) is certainly not larger than that caused 
by the use of mean opacities in the source moments (see below) and, hence, is adequate for the present 
frequency-integrated calculations. 

In this paper we will use very simple input physics for the coupled gas equations of state and the source 
moments. For our simplest models we will adopt a totally ionized hydrogen gas with constant opacity 
(independent of p and T). From the numerical calculations we have found that, during computation, the 
maximum gas temperature is reached near the horizon and is not larger than ~ 5 x 10 8 A'. In our more 
detailed models we will consider a hydrogen gas with variable degree of ionization, computed using the 
Saha equation. In a forthcoming paper we will treat a more realistic gas composition. If the electrons are 
nonrelativistic (T < 10 9 K), we have 

p=[l + x{T )]PML, (14) 

e = p{l + | [1 + x(T)] ^ - [1 - x{T)] — \ > (15) 

where T is the gas temperature, Ep, = 13.6 eV is the hydrogen ionization potential and x(T) is the degree 
of ionization, given by 

X ^ = TT5- (16) 

We assume that all atoms are in their ground state. The function 5 is the positive root of the Saha quadratic 
equation (Rybicki & Lightman 1979) 

^ 2 _ m H ( 2-Km e k B T \ 3/2 c - Ell /k B T ( 17 ) 
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Equation ( |17| ) holds in the "one level atom" limit for hydrogen, i.e. when the temperature is fairly low 
compared to Ep. Although this is not going to be true for our calculations throughout all the grid, however 
it is still true that equation ( |17j ) will give a large ionization at high temperature. The third term inside 
the curly brackets in equation (|l5|) accounts for the electrostatic potential energy of bound electrons in the 
neutral hydrogen atoms. For the opacities, we adopt the interpolation formula by Christy (1966) which, for 
the range of temperatures and densities considered, fits the hydrogen Rosseland mean opacity kp to better 
than 30%: 

1 i/ 2 /2xl0 6 



fc^ = Pe ^4.85x 10- 13 — + T^-jj— +2.13? 
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(18) 



where p e is the electron pressure and T4 = T/IO^K. The Planck mean opacity kp has been calculated 
approximately by subtracting the Thomson scattering contribution (the first term in the curly brackets 
in equation |jl8|| ) from the expression for kp. This choice follows from the fact that a very accurate 
evaluation of kp is not critically important since the transition of the ejecta from optically thick to thin 
is governed entirely by the recombination process. Furthermore, no detailed analytic interpolation of the 
Planck mean opacity in the range of temperature and density of interest here can be found in literature and 
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numerical interpolation from accurate tabulated opacities would be too expensive in the present preliminary 
calculation. Using the Kirchhoff law, the first two radial moments of the source function so and s\ can be 
written 

s = p(k P B - fc W(i) , (19) 
81 = -pkxwi , (20) 

where B — a^T 4 (ap is the blackbody radiation constant) and fco & n d k\ are the absorption and flux mean 
opacities. Since in the present frequency-integrated calculation the actual spectral distribution of wq and 
w\ is not known, we use kp and kp in place of kg and fci, respectively. Finally, the optical depth r at radius 
r is defined by 

/•OO 

t= k R pdr. (21) 



3. Numerical method and boundary conditions 

The numerical code used to solve the equations of radiation hydrodynamics and the radiation moment 
equations presented in the previous section is based on a Lagrangian finite difference scheme with a standard 
Lagrangian organization of the grid (see the Appendix; see also Zampieri, Miller & Turolla 1996). The 
spatial and time centering of the variables ensures second order accuracy both in space and time. As far as 
the spatial centering is concerned, p, B = apT 4 , wo and a are evaluated at mid-zones, while r, M, u and 
wi are evaluated at zone boundaries. To have second-order accuracy in time, u and w\ are both evaluated 
at an intermediate time level (time-shifted). Because of the way in which variables are centered, the code 
is semi-implicit. This feature helps preserve stability (see e.g. Mihalas & Mihalas 1984; Shapiro 1996). The 
full set of equations in finite difference form is written in the Appendix. 

The time step is controlled by the Courant condition for stability and by additional constraints on the 
fractional variation of the variables (< 10%). The Courant condition reads 

At 

aTv c — < 1 , (22) 
Ar 

where v c = (/ + 1/3) 1 ' 2 is the characteristic speed of radiation in the Lagrangian frame (the fastest 
characteristic speed on the grid in the present problem). The grid in Lagrangian mass \x is usually divided 
into jmax = 200 zones. A constant fractional increment in grid spacing between successive zones a is used. 
It is calculated from the equation 

Jmo»-l 

Vjmax = fJ'jmin + A Vj + l/2 > ( 23 ) 

where fJ>j min is the inner boundary, fXj max is the outer boundary (fixed for the conservation of the total 
envelope mass) and A/ij +1 / 2 = Mj+i ~ Mj = aA Pj-i/2- At the beginning of the calculation fJij mi „ — 0. The 
mass contained within the first shell A/ix/2 is fixed by the requirement that the radial spacing between the 
first two shells is 30% of rj min . Equation (|2^) can be written as a polynomial of order j max in a and has 
been solved using the Newton-Raphson method. Whenever the inner edge of the innermost zone crosses 
the inner boundary in radius rj n , it is removed from the calculation and a regridding of all the variables is 
performed. During computation, the actual inner boundary in mass changes according to: (J-j min = (J>j min +u 
where I is the number of zones that have crossed 7"j n . The fractional increment a is computed solving 
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equation ( p3| ) with the new value of fJ>j min ■ All the variables are then interpolated on the new grid. The 
regridding procedure makes use of a local cubic interpolation that, for uniform spacing, reduces to a local 
fourth-order Lagrangian interpolation. Performing systematic regriddings is necessary both to minimize 
the number of grid points and to preserve the spatial resolution with time. 

Once the finite difference representation has been introduced, equations (^), (|J), (||), (Q) and (||) can be 
solved explicitly (algebraically) for u, p, a, M and r, respectively (see Appendix). Where necessary, linear 
interpolation and extrapolation in time were used to obtain the values of quantities at suitable time levels. 
Using equation (20) with k\ — Ur, equation ( fill) is rewritten in the following form 
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and solved (algebraically) for wi at level n + 1/2. The quantity W2 is expressed using the closure relation 
u>2 = fu>a with / being defined as in equation (|l2|). Again, where necessary, variables on the right hand 
side have been interpolated or extrapolated at the correct time level (see Appendix). The energy equation 
(0) and the 0-th moment equation (M) form a strongly coupled, nonlinear system of equations. Using 



equations (19) with k Q — hp, they can be cast into the form 
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«to + -^(wiaV), /1 = 1 (26) 



where e = (e — p)/ p is the internal energy per unit mass. They have been solved for B and wo at the new 
time level n + 1 using the Newton-Raphson method for nonlinear systems of equations. 

As far as the boundary conditions are concerned, the time slice at constant t is a characteristic direction 
for equations (||) and (0). At the outer boundary we put 



1 



(27) 



Equation ( p7|) corresponds to synchronizing the coordinate time with the proper time of a comoving 
observer at the outer edge of the grid. The inner boundary is chosen away from the origin to avoid the 
black hole singularity. There we set 



M = Mi 



o 



/- 1 = Hmin > 



(28) 



where Mq is the effective mass contained within the inner boundary. At t = this is taken to be equal to 
the initial black hole mass Mt,h- I n the limit appropriate here where the infalling fluid at Pj min is highly 
supersonic with negligible thermal and radiation back pressure, the gas behaves like dust. In this case, as 
time increases, we know that the effective mass M increases at the rate at which baryons (i.e. rest mass) 
cross into rj„. Once they cross, they are removed from the calculation. The boundary conditions for the 
continuity and Euler equations are 



u.t 
u.t 



aM 











(29) 
(30) 



Equation ( p^ ) corresponds to assuming that pressure gradients and radiative forces can be neglected at the 
inner boundary and that the free fall is achieved near to /i = p.j min ■ This turns out to be correct if the inner 
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boundary radius rj„ is sufficiently smaller than the accretion radius (r a = GMbh/c^). During the initial 
evolutionary phases, in some of our computed models radiative forces are very important and affect the flow 
dynamics. However, in the innermost part of the accreting region the luminosity is always sub-Eddington 
(see Section^). Thus, if the inner boundary is sufficiently far in, at /i = fij min the gas is moving in near 
free-fall. The condition at fi = [ij max follows from the requirement of free expansion. It is correct as long 
as the motion of the gas near the outer boundary is supersonic and both radiative and gravitational forces 
have little influence on the outflow (see e.g. Arnett 1980). The inner boundary condition at /i = [ij min for 
wi is a five-point Lagrangian extrapolation in r along the time slice at constant t. To impose the boundary 



condition at the outer edge of the grid we apply equation (|24| ) across the half-zone from j max — 1/2 to 
jmax, substituting (wi)j max = g(wo)j mam for (wo)j max in the gradients with respect to /z (see e.g. Mihalas & 
Mihalas 1984; Shapiro 1996). The closure boundary factor g is given by equation (|l3|). 

As already mentioned, the evolutionary times for the present problem are very long in comparison with 
the timescales involved in the calculation. For this reason we developed a multiple timestep procedure to 
speed up the calculation, discussed in the next section. 



4. Multiple timestep procedure 

A situation often encountered in the numerical solution of time-dependent problems is the need to 
integrate the equations over a very wide range of length scales and corresponding timescales. Frequently, 
the relevant timescales vary monotonically along the grid. In these situations, integrating every cycle over 
all of the mesh with the smallest timestep is an unwanted waste of computational time (in particular if 
the number of grid points is large and/or the evolutionary times are very long), since the evolution is 
unnecessarily slow in some portions of the grid. We have developed an approach, that seems particularly 
suitable for the present problem. The basic idea of our multiple timestep procedure (MTP) is to decompose 
the grid into subgrids, each one evolving separately according to its own timestep. The subgrids with 
larger timesteps are evolved with fewer steps than those with smaller ones. The major challenge here is 
establishing the communication between neighboring subgrids. Additional complications arise if the grid is 
staggered in time, so that some variables are evaluated at the full time level (time centered) whereas other 
variables are evaluated at half time level (time shifted). Finally, if the code performs systematic regriddings 
over all the grid, the evolution has to be re-synchronized before performing any interpolation. Our code 
handles all of these difficulties. 



4.1. Communication at the boundaries between neighboring subgrids 

In the following we assume that the integration domain has been divided into N subgrids. We assume 
also that the innermost zone (k = 1) has the smallest timestep and that the sequence of subgrids has 
timesteps monotonically increasing outwards. This assumption is certainly correct in the present case. In 
deciding the order in which to evolve the various zones, we must consider the subgrid boundaries. Subgrids 
have boundaries which are in the interior of the integration domain and hence they will need boundary 
values to close the PDE's. In the Adaptive Mesh Refinement method a coarse grid is advanced (in time) 
first and then all the subgrids are integrated to the same time level. Boundary values are often computed 
by interpolation from the coarse grid (Berger & Oliger 1984). By contrast, we do not evolve a coarse grid 
and hence we cannot interpolate variables at the intermediate boundaries. Instead, we evolve the various 
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subgrids (zones) in order of increasing timesteps (i.e. in the sequence k = 1, k = 2,..., k — N). In this way, 
all of the subgrids, apart from the last one (k = N), need values of same variables, extrapolated from earlier 
times, at the outer subgrid boundary. On the other hand, boundary values at each subgrid inner boundary 
can be computed by interpolation in time from values previously calculated in the adjacent subgrid. Linear 
extrapolation and interpolation in time has been used. 



4.2. Subgrid evolutionary algorithm 

The subgrid evolutionary algorithm establishes the order and the conditions under which the various 
subgrids have to be advanced in time. The main criterion is to evolve the subgrids in order of increasing 
timestep. Initially, the solution is evolved for a given number of cycles n g over all of the grid with the 
smallest timcstcp. When n > n g , the subgrid evolutionary algorithm is activated. The evolutionary time 
t^ k and the timestep At^ fc+1 ^ 2 of all the subgrids are stored in a vector that is periodically updated. The 
choice of the subgrid evolved at each cycle and its timestep are decided according to the following steps: [1] 
The evolution always starts from the innermost [k = 1) subgrid, so that i™ 1 > t^ 2 ■ The innermost subgrid 
is always evolved according to the Courant timestep, which takes its minimum value in this zone. [2] The 
innermost subgrid is advanced in time until the difference between its coordinate time and the one of the 
neighboring subgrid does not exceed a given threshold: 

A/ «2 + l/2 

t? - tT < (31) 

where 5 > 1 is a fixed number. Inequality ( |3l| ) tells the code if the difference in the evolutionary times 
between subgrid k = 1 and k — 2 is less than S^ 1 times the timestep of subgrid k = 2. [3] If the difference 
between the coordinate times exceeds the threshold, subgrid k — 2 is advanced to the time level of the 
first subgrid so that t£ 2 = t" 1 . At this point the time centered quantities of subgrids k = 1 and k = 2 are 
synchronized (but not the time shifted). [4] Then, the code checks whether the inequality corresponding to 
( |3l| ) for subgrids k — 2 and k = 3 is satisfied: t% 2 — t^ 3 < At 7 ^ i+1 ^ 2 /S. If the separation in time is below 
the threshold, the evolution continues from the first subgrid. [5] If it is not, subgrid k = 3 is advanced 
up to the time level of subgrid k = 2 (which, at this stage, is also equal to t" 1 ). Time centered variables 
are then synchronized (whereas time shifted quantities are not evaluated at the same time level). This 
procedure is repeated for all the subgrids. Note that the factor S has a very important role. The smaller S 
is the faster the code will evolve, since the A;-th subgrid will be advanced fewer times with respect to the 
(k — l)-th subgrid. However, too small a value of 6 may cause accuracy and stability problems because 
of the extrapolation of the variables at the outer subgrid boundary. By experimenting, we found that the 
value 5 ~ 5 turns out to be a good compromise between computational speed and accuracy. 



4.3. Synchronization at regridding 

The innermost subgrid is regridded whenever the inner boundary of the innermost mass shell rj min 
crosses radius rj„. Synchronization at regridding represents a major problem since not all variables on the 
grid are time-centered. We solve the problem in the following way. Suppose that we are in the need of 
performing a regridding and let N be the number of grids which are advanced up to the same time level 
(see previous subsection) . All the time-centered quantities of the subgrids interior to (and including) N will 
be synchronized. As described in Section 13, we need to interpolate in radius all the variables (time — shifted 
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and time-centered). This means that also the time-shifted variables have to be synchronized. This is 
accomplished evolving all the k < N subgrids together for two cycles, i.e. forcing the timestep of all the 
subgrids to be equal to the Courant timestep of the innermost region. In this way, both time centered 
and time shifted quantities in all the subgrids with k < N are synchronized and a regridding of all of the 
variables can be performed. If N — N, all the mesh is regridded. 

4.4. Performance 

We have tested the MTP procedure in different situations and using the test problems which will be 
described in the next sections. We have varied different parameters (S, N, j max , the number of points of 
each subgrid, etc.) and the procedure turns out to be stable and accurate. The gain in computational time 
for a run with a grid of 200 zones in a domain spanning 4 decades in radius is about a factor 7-8 using 
4-5 subgrids. We found that, beyond a certain number of subgrids (which depends on the specific problem 
and integration domain), the gain in computational time depends mainly on the number of points and 
the position of the boundary of each subgrid and not very much on N. A typical run including hydrogen 
opacity and recombination takes about 10 days on a Digital ALPHA station 200 4/233. 

The main limit of the MTP procedure is certainly related with the extrapolation to get boundary 
values at the subgrid boundaries. This may cause problems following the evolution of shocks or any type of 
sharp feature possibly present in the solution. Problems may also originate if the flow enters a regime in 
which there are very delicate balances of forces, as for example when the luminosity is very close to the 
Eddington limit (within 1%). In such cases, it is necessary to restore the one timestep integration for the 
time necessary to follow the development of sharp features or the evolution of sensitive regimes. However, 
if the prominent feature is not crossing the boundary between two subgrids or if the delicate regime is 
contained within one subgrid, the evolution can proceed reliably using the MTP procedure. 

5. Initial conditions and relevant timescales 

In a supernova explosion, fallback of material may occur when a reverse shock forms at the interface 
between the mantle and the hydrogen envelope of the progenitor star (see e.g. Woosley & Weaver 1995). 
After the shock has propagated toward the center, the mantle is set into homologous expansion and matter 
bound to the central remnant is accreted. A complete analysis of the fallback process should comprise the 
simultaneous study of the expanding hydrogen envelope and of the mantle, a region where large composition 
gradients develop. Due to the complexity of the process here we explore a simplified model that follows 
the radiative and dynamical evolution of an expanding "cloud" of fixed composition. All of the emitted 
energy comes from the release of heat residing in the gas originally or is generated by comprcssional heating 
in the course of accretion onto the black hole. Radioactive energy sources have not been included in our 
calculations. Thus, our model will not give an accurate representation of the light curve of a supernova 
which is provided, at and after the peak, by radioactivity. Our main purpose is simply to calculate the 
contribution coming from accretion onto the central black hole and give quantitative estimates of its 
importance. 

At the onset of evolution the cloud has homogeneous density po and is set into homologous expansion 
with a velocity profile u = r/to, where to denotes the expansion time scale. The maximum velocity of the 
ejecta is Vq = u(r out ) = r out /t 0} where r out is the outer radial boundary of the expanding cloud. The 
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temperature profile is taken to be equal to the "radiative zero solution" of Arnett (1980) 

r • t m i/4 
sm(7nE) 



T(r,0)=T o 



where T is the initial temperature at the inner boundary and x = r/r a 



(32) 



Four parameters specify uniquely the dynamical and thermal state of the cloud, for a fixed composition 
the total mass M c i ou d, the radius r out , the sound speed at the inner boundary c Si o, and the ratio k of the 
accretion timescale t a ,o (defined below) to the expansion time to- Five relevant timescales are involved in 
this phenomenon. Using the electron scattering opacity (for completely ionized hydrogen) as reference 
value, k es — 0.4 cm 2 g _1 , they are defined as follows: 
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(35) 
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t a fi an d to are the initial accretion and expansion timescales. The diffusion timescale, tdiff, represents the 
time for the radiation to diffuse out of an homogeneous expanding cloud. It can be estimated considering 
tdiff ~ T es (tdif f)r ou t(tdif f) I '3c, where 3c/T es (tdiff) is the photon diffusion velocity and T es is the electron 
scattering optical depth. The timescale t tra ns represents the time at which the cloud becomes transparent 
because of the decrease in density caused by the expansion; it is the time at which the photon mean 
free path \(t tr ans) = l/k es p{ttrans) is of the order of the cloud radius r out (i 4r£ms ). For free expansion, 
the density decreases as p oc t~ 3 . The recombination timescale t rec is estimated calculating the time at 
which the gas temperature drops below the recombination temperature T reci as a consequence of adiabatic 
expansion; for a T = 4/3 polytrope, T oc p 1 / 3 oc 

We consider a series of models selected on the basis of a predetermined hierarchy of timescales which 
will guide our understanding (for the parameters and related timescales see Table 0). Models 0, I, II and 
III deal with a cloud of fully ionized hydrogen for which the recombination time t rec — > oo. Models IVa and 
IVb allow for recombination: they have T rec ~ 10 4 i^ and t rec < tdiff. In all models, the mass of the central 
black hole is taken to be M^h = 1.5M Q . 

The first set of models (0 and I) have been introduced mainly to test our general relativistic radiation 
hydrodynamic code against previous numerical and analytical results. The choice of the parameters for 
the other models follows from the requirement that the hierarchy of relevant timescales reflects that for a 
realistic supernova model. For example, for the inferred parameters of the hydrogen envelope of SN1987A, 
we have: t a ,o <C to -C t rec < tdiff <C t tra ns (see Table ||). The second set of runs (models II and III) do not 
include recombination and for them t a ,o to *C tdiff ttrans <C t rec . In practice, in these models, the 
photons have time to diffuse outwards before recombination takes place. Although they do not reflect the 
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true physics, these models are important to understand the role played by the central accreting black hole 
in shaping the light curve (in particular at late time). Since they have mainly a "pedagogical" purpose, 
we decided to run these models with a relatively small cloud mass (M c i ouc i = 10 _2 Mq) and in turn small 
accreted mass. In this way and t trans are decreased and, consequently, the computational time is 
significantly reduced. Note that, because tdiff and ttrans have the same dependence on the cloud mass, 
decreasing M c i oua - does not alter the hierarchy of timescales as long as to < tdiff' It should be noted also 
that in varying k (i.e. the expansion velocity Vq) the hierarchy of "radiative timescales" remains unchanged. 
In fact, from equations (|35|) and (E3) 



t d if L = ±(Vo} 1 ' 2 



ttrans v3 \ C 

In examining models I, II and III, we will assume that the opacities are given by 



k P = 10" 2 A;e 6 
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(39) 



k R = k es = 0.4 cm 2 g- 1 . (40) 

The expression for k has a maximum at r = rj min and then goes rapidly to a constant value of 10 _2 fc e s- 
This expression has been used for convenience to keep the inner core in LTE as long as possible in order to 
avoid the decrease in time-step caused by the inward motion of the inner boundary. 

The last set of runs (models IVa and IVb) is a step forward to the study of more realistic models 
with t a ,o <C to <C t rec < tdiff <C ttrans- Table H contains the relevant parameters for SN1987A describing 
the post-shock structure of hydrogen envelope and mantle (after the reverse shock has reached the center) 
as deduced from Arnett (1996) and Chevalier (1989), respectively. As can be seen by comparing Tables 
|l| and |[ model IVb has quite realistic initial parameters, although they do not match exactly the ones 
characterizing the hydrogen envelope in SN1987A. The calculation of a SN1987A-like model would be at 
least 5-6 times longer computationally. Because of the preliminary nature of the present investigation, we 
have not attempted such a calculation. 



6. Numerical results 

In this section we present results of the numerical calculation of the dynamical and thermal structure 
of the expanding cloud, starting from the initial conditions outlined in the previous section. All the light 
curves presented in the figures have been computed using u>i calculated at the last grid point [ij mam (which 
corresponds to the outermost radius r out ). It is 

L{r out ) = Airr 2 out wi (r out )c . (41) 

In equation (|4l|), L(r ou t) is the luminosity measured in the comoving frame. To obtain the luminosity 
seen by a distant Eulerian observer L 00l one must transform to the Eulerian frame. However, at the outer 
boundary r out , general relativistic effects can be neglected and the flow velocity Vo never exceeds 0.1c. 
Then, to calculate we can safely ignore relativistic corrections and write 



Loo — L(r out ) 



(42) 
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In the following we will always refer to the luminosity and accretion rate in units of the corresponding 
Eddington values for Thomson scattering in a completely ionized hydrogen gas. It is 

I = -A- (43) 

i-'Edd 

m = (44) 

Msdd 

where L Edd = Ai:GM bh c/k es = 1.3 x 10 38 (M 6/t /M s ) erg s"\ M Edd = L Edd /c 2 = 2.3 x 10" 9 (M bh /M Q )M Q 
yr -1 and the local accretion rate Mq has been defined by 

AIo = 4-rrr 2 puc . (45) 



6.1. Model 0: polytropic accretion comparison with CSW 

These models have been run to test the code and to compare our numerical results with those obtained 
by CSW. Two set of calculations were performed: accretion of a polytropic T = 4/3 gas with no radiation 
(see also Zampieri 1997) and accretion of a radiation dominated, completely ionized hydrogen gas with 
constant opacity. In the first set of runs, we considered two models with k = 0.01 and k = 1 and with the 
same value of the initial sound velocity (c Si o = 10 8 cm s _1 ). The agreement with the results of CSW is 
good: the density profile and the accretion rate as a function of time differ at most by ~ 5% percent. As 
already mentioned, we have tested these models changing a number of parameters such as the number of 
points, the position of the inner boundary, the number of subgrids, the position of the subgrid boundaries 
and S. We found that the numerical results agree to within a few percent. The second set of runs includes 
radiation. The initial parameters are the same of the polytropic models. The absorption opacity was taken 
to be sufficiently high that, during dynamical evolution, matter and radiation were always in LTE over 
most of the integration domain. Furthermore, the evolutionary time is smaller than the diffusion timescale, 
so that the evolution is quasi adiabatic. In these conditions we expect that the solution should agree with 
the one obtained in the adiabatic case. The comparison of the accretion rate versus time with the results for 
a r = 4/3 polytrope shows agreement at the 1-2% level. We have also tested that, rescaling the flux mean 
opacity, the flow dynamics remains unchanged. This is the expected behavior for a radiation dominated 
optically thick gas, since the radiative force kiwi oc (l/3&p)(u>o).p and the radiation energy density gradient 
is basically insensitive to an overall rescaling of fci. 



6.2. Model I: comparison with Arnett's analytic solution 

If the initial accretion timescale t a ,o is longer than the diffusion timescale tdiff, the thermal and 
radiative structure of the envelope resembles that of a hot cloud that is cooling because of expansion and of 
emission of radiation. Figure [l] shows the light curve for model I for which t a fl ~ tdiff- The filled squares 
are the computed points and the solid line gives the analytic solution for a pure radiation pressure, constant 
density (and opacity) gas given by Arnett (1980; "radiative zero solution") 

I = / oe -[*/td,//,o+t 2 /(" i t2,//)] , (46) 
where tdifffi = 3k es por 2 ut /ir 2 c, ai — 18/tt 2 and Iq is the luminosity scale 

— ^f- (47) 
3 r s p Q c/ 
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In equation (^7j) r s is the Schwarzschild radius. Since idi//,o ^> tdiff, at late times the surface luminosity 
falls as a gaussian. As can be seen from the figure, the agreement between the numerical and analytical 
results is excellent. For t < t^ift the luminosity is constant because the decrease in the radiation energy 
density gradient is exactly compensated by the increase in both the photon mean free path and the 
photospheric radius. When t becomes larger than tdiff, the internal energy is being exhausted and the 
luminosity falls off as a gaussian. 

6.3. Model II: photon diffusion in a high temperature, expanding cloud with a central black 

hole 

For the typical parameters of a supernova with a black hole mass of ~ 1.5M Q , the accretion timescale 
is usually much smaller than the diffusion timescale. The presence of the central black hole is "felt" by the 
expanding cloud. We have computed a model with M c \ ou d = 10~ 2 M@, corresponding to tdi///i o ,0 ~ 5. 
The results are shown in Figures In Figure [|, for t < tdiff the luminosity is roughly constant and is 
still due to the emission of the cloud internal energy. The differences with respect to the analytic solution 
(equation [^6|) can be explained in terms of the different hierarchy of timescales (t Qj o < tdiff) and of the 
increasing importance of the gas pressure. In fact, for this model the initial value of the radiation to gas 
pressure ratio is not very large (in particular near to the outer boundary where the temperature drops 
rapidly). Then, equation (|4^) is no longer a good approximation to the numerical solution. At t ~ tdiff, 
the initial exponential decline in the light curve turns into a steep power law. This is the first signature 
of the presence of the central black hole. In fact, if no black hole were present, one would expect the light 
curve to fall off as a gaussian with time. The radiative energy emerging at this stage has been produced by 
accretion during the early phases of the evolution and transported outward diffusively. At ro = 0.3r aj o the 
maximum value of the diffusive accretion luminosity I ~ 0.5 is reached at t ~ i aj o (see Figure |^a) slightly 
after the peak in the accretion rate (Figure ||b). After traveling outwards through the opaque cloud, this 
accretion energy emerges on a diffusion timescale: this corresponds to the point, along the light curve, at 
which the analytic calculation fails to reproduce the actual behavior. Later, the cloud becomes transparent 
(t ~ ttrans), and the central accreting black hole "becomes visible". The steep decline of the light curve 
ceases and the emergent luminosity is entirely produced by accretion of the gravitationally bound gas. The 
internal accreting region is still optically thick. We note that, at times before tdiff, the heat content of 
the gas is dominated by radiation, but at times after tdiff it becomes dominated by the electrons, protons 
and hydrogen atoms. At this stage, the timestep of the numerical integration becomes exceedingly small 
because the inner boundary has to be moved inward to recover the region where the radiation field is in 
LTE. In practice, this prevents us from studying numerically the very late-time evolution of the model. 
However, as we shall see later, the light curve can be equally well computed using a sequence of stationary 
models with progressively decreasing density at infinity (see next section). 

Figure ^ summarizes the properties of the flow at different times along dynamical evolution. As shown 
in Figure |Ja, at the onset of expansion the gas is outflowing with u cx r. Very soon the inner low velocity 
shells start to accrete and free fall (u cx r -1 / 2 ) is rapidly achieved. Because the dynamical timescale in 
the inner part of the flow is much smaller than that in the outer envelope, in the inner accreting region 
the density profile adjusts very quickly to the p cx r~ 3 / 2 behavior (Figure [|b), as in stationary flows. In 
the outermost regions p remains nearly independent of radius and decreases with time (~ t~ 3 ) because 
of expansion (CSW). The evolution of the gas temperature is closely coupled with the properties of 
the radiation field (Figures ^c and^jd). In the inner accreting region, the temperature increases inward 
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(T oc r~ 5 / 8 ; see e.g. Blondin 1986 equation [17]) because of the compressional work done by gravitational 
forces. In the outer expanding region, at early times T remains roughly constant with radius, except close to 
the outer boundary, where it drops significantly to transport photons outwards diffusively. In this respect, 
the temperature profile reflects the behavior of the solution of Arnett (1980). As shown in Figure ||c, during 
early evolution when gas and radiation are in LTE and t < tdiff, T drops adiabatically with time (T cx t^ 1 ), 
as expected for a freely expanding radiation dominated gas. In Figure |]d, we have plotted the luminosity 
profile at different times along dynamical evolution. At t ~ 3i ai o, I has two maxima: the first represents 
the accretion luminosity produced by the central black hole, whereas the second is due to the emission of 
the internal energy stored in the cloud. After the initial transient phase, the accretion luminosity decreases 
steadily with time. The diffusive luminosity of the cloud stays almost constant up to t ~ tdiff and then 
falls off as a gaussian. At t ~ 500i aj o the cloud becomes transparent to the central accreting region (see the 
curve at t = 1034i Qj0 in Figure §d). 



6.4. Model III: photon diffusion in a high temperature, slowly expanding cloud with a 

central black hole 

As we already noted, the hierarchy of "radiative" timcscales remains unchanged as we vary the velocity 
of the ejecta Vq. So the light curve of model III (same parameters of model II but k — 0.1; see Table |l|) 
shows the same basic features of model II: initial plateau due to the emission of internal energy, a steep 
exponential decline turning into a power law at t ~ tdiff because of the emergence of the initial transient 
accretion energy, cessation of the rapid decline due to the clearing of the cloud, revealing the presence of 
the central black hole (see Figure ||). However, some differences can be noted. First, comparing Figure 
|5] with Figure 2, it appears that the power law decline after t = tdiff has a different slope. This can be 
explained in terms of the different expansion timescale. Since model III has an outward expansion velocity 
10 times smaller than model II, less work is being done by the gas and the internal energy suffers less 
adiabatic degradation. Also the energy released by accretion undergoes less degradation and can then 
produce a flatter tail in the light curve after t ~ tdiff. Second, because of the lower expansion velocity of 
model III, for a certain time the outer part of the expanding cloud remains opaque whereas the central 
part is optically thin. During this phase, the accretion luminosity produced in the inner region is partly 
absorbed by the outer optically thick region. This is the reason for the decrease in luminosity close to the 
outer boundary at t ~ 2400i a .o in Figure ^. 



6.5. Model IV: clearing from expansion and recombination 

As illustrated in Figure ^c, by the time that the cloud becomes optically thin because of expansion the 
temperature has already dropped well below 1Q 4 K, the temperature at which hydrogen recombines. In fact, 
in a realistic model of a supernova explosion, the clearing of the hydrogen envelope takes place because 
of recombination. For this reason we consider a pure hydrogen gas with variable degree of ionization in 
model IVa and IVb. The ionization coefficient is computed from the Saha equation. As can be seen looking 
at the temperature and luminosity profiles of model IVa (Figures |a and |b), for t < t rec the evolution is 
very similar to that of model II. At t ~ t rec hydrogen starts to recombine. The recombination process is 
very fast and generates a recombination wave that propagates rapidly through the envelope (Figure |^a) . 
Outside the recombination front the gas is optically thin, whereas inside the mean free path is so small that 
the gas is everywhere opaque. The recombination front almost coincides with the photosphere. Note the 
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sudden rise in luminosity in Figure |8|b after the recombination front has propagated inward. As noted by 
Woosley (1988), during this phase, the radiation does not have time to diffuse through the photosphere; 
the photosphere instead moves to the radiation, releasing all of the internal energy. The huge amount 
of internal energy released gives rise to a big bump in the light curve at t ~ t rec (see Figures ^ and ^a) . 
After the maximum, the luminosity falls off abruptly (see also Figure ||b). In contrast to models II and 
III, no power law decay can be observed after maximum. The reason is that the fast inward motion of 
the recombination front has liberated all of the internal energy before it could diffuse outwards. Also the 
initial transient accretion energy, stored as internal energy in the gas flow, has been completely liberated 
and contributes to the big recombination bump in the light curve. During this phase our computed light 
curves are similar to those calculated by Woosley (1988) for stars with low envelope mass. In this respect 
they resemble the light curve of a Type II-Type lb supernova. At late times, after that the recombination 
front has propagated down to the inner accreting region, the light curve is entirely powered by accretion 
onto the central black hole and the luminosity decreases as a power law with time (see next section). In 
model IVb, the total mass accreted onto the central black hole is ~ 2 x 1O^ 5 M . 

In Figure ||c we follow the evolution of the recombination front as a function of time (for model IVa) . 
The front clearly stalls as soon as it approaches the innermost region where compressional heating due 
to accretion overcomes radiative cooling. The density and velocity profiles are largely unaffected by the 
propagation of the recombination front, which instead influences the light cure and thermal state of the 
cloud. In Figures ^|b we plot the effective temperature, T e ff 7 and the gas temperature at the photospheric 
radius, T p h, as a function of time for model IVb: 

Teff = ^0 T ph = T(r ph ) , (48) 

where r p h is the photospheric radius (i.e. the radius at which the effective optical depth is equal to unity). 
As can be seen from Figure [)]b, during the phase of "adiabatic" expansion, the effective temperature 
decreases with time to ~ 4000/4'. When recombination takes place, the photosphere starts to recede through 
the envelope and T e // rises to ~ 1Q A K. From this moment on the photosphere is located in the inner 
accreting region where the compressional heating is balanced by radiative recombination. Furthermore, for 
t > 10 3 £ a .o the gas photospheric temperature stays very close to the effective temperature, which suggests 
that the emitted continuum spectrum should be roughly a blackbody. When all of the outer hydrogen 
envelope is recombined and all of its internal energy has been emitted, the fall in luminosity ceases and 
the light curve starts to be dominated by the accretion luminosity of the gravitationally bound shells. At 
this point the timestep of the numerical integration becomes very small. However, as we will discuss in the 
next section, the asymptotic curve I = l(m) can be estimated using a sequence of stationary models with 
progressively decreasing density at infinity. 

Consider the early behavior of model IVb (see Figures |To|), which has the largest initial cloud mass 
M r c ioud- The large amount of mass that falls back within the first few hours gives rise to an initial (diffusion) 
luminosity transient that reaches the Eddington limit (see Figure [To|b) . At this point the flow starts to be 
"self-regulated". As the accretion rate increases, it drives the luminosity above the Eddington limit. Then, 
the radiative force eventually pushes outwards the accreting gas (note the sign inversion of the velocity 
of the intermediate shells at t ~ 2t 0) o in Figure |Io|a). The flow readjusts itself in such a way to maintain 
the luminosity close to the Eddington limit. This delicate balance between radiative and gravitational 
forces governs the flow dynamics and allows the gas to radiate very close to the Eddington limit during 
the first evolutionary phase. This result is in agreement with the findings of CSW. In their T = 4/3 
polytropic calculation, they have shown that in flows with the pressure pushes initially bound shells 
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outwards, preventing them from accreting. Since an optically thick, radiation dominated flow is dynamically 
equivalent to a T — 4/3 polytrope, the gas pressure gradient in the CSW models acts much in the same way 
as the radiative force in our models. 



7. Late time light curve 

CSW showed that in polytropic flows with k = i a ,oAo > 1, pressure forces act as a small perturbation 
and the accreting gas behaves as a nearly collisionless fluid. In this case the total mass accreted is (CSW, 
equation [17]) 

&7T 

M dust ~ — GM bhPo tl , (49) 

and the late-time evolution of the accretion rate (in units of the Eddington rate) is equal to (CSW, equation 
[29]) 

4vr 2 / 3 / t \ " 5/3 

in{t)— q Potock es i—j . (50) 

In the opposite regime, k = i aj o/io < L pressure forces significantly affect the motion. After a time 
t ~ 10t ai o the flow evolves along a sequence of quasi-stationary, Bondi-like states, with a time dependent 
m determined by the slowly varying density at large distances (CSW, equation [45]) 

m(t) ~ \ r GM bh ck es -^P- , (51) 

where c Sitx) = (TKp^ 1 ) 1 / 2 is the sound speed (K and T are the polytropic constant and the polytropic 
index) and A r = 0.25[2/(5 - 3r)]( 5 - 3r )/2(r-i)_ For homologous expansion Poo = p (l + kt/t a _ )- 3 . The 
persistence of a Bondi-like behavior is however limited to an interval of time between 10t a fi < t < t tr . 
At late times t > t tT1 the flow enters the dust regime and the accretion rate falls like m cx t~ 5 / 3 . The 
transition time tt r is estimated to be ~ 9fc~ 3 i aj o/2, for T — 4/3, and is the time when the instantaneous 
Bondi accretion radius reaches a position where the fluid is unbound and flying away supersonically. For 
a flow with k < 1, the mass that has fallen back can be estimated multiplying the Bondi accretion rate 
MBondi — 47rAr/Oo r a oAa,o by the accretion timescale t Qi o, which gives 

M Bondi ~ 4^A r p [GM & b ^ ■ (52) 

c s,0 

As can be seen from the previous equation, the total accreted mass is very sensitive to the actual value 
of c St o. MBondi is only a tiny fraction 3Ar£; 2 /2 ~ t\ o/tg of the mass accreted if the flow was collisionless 
(see equation f49f| ). In fact, for fixed po and to, the most favorable situation (when the accreted mass is 
maximum) occurs if the sound velocity is so low that pressure forces can be neglected, i.e. if the flow is 
collisionless. In this regime, to becomes the relevant timescale for accretion and the total accreted mass is 

~ Mdust- 

Models IVa and IVb have k = 1 and k — 0.1 respectively. In IVa, the fluid is nearly collisionless 
and, after time t ~ 10£ Qj o, we find that fn (plotted in Figure |li"| ) is well approximated by equation (|5C|). 
In IVb, the expansion timescale is instead larger than the accretion time and, at intermediate times 
(t a ,o < t < 4.5 x 10 3 t Qi o), the flow is expected to evolve along a sequence of quasi-steady states characterized 
by T = 4/3 and m cx t~ z l 2 (equation J5l|). As shown in Figure O, the accretion rate decays approximately 
as i -3 / 2 , as long as t <; 300t a .o- We then observe a slight bending toward higher rates. Later on, the 
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accretion rate changes slope again and rh starts decaying as i~ 5 / 3 : thereafter, the flow behaves as dust. 
Note that, due to the small envelope mass (for IVa) and small k (for IVb), the total mass accreted is tiny 
in these two models. 



In Figure 12 we have plotted rh as a function of radius at t — 6 x 10 3 i Qi o for model IVb. As can be 
seen, in the inner region the accretion rate is nearly independent of radius. This means that the evolution 
follows a sequence of quasi-stationary states. Indeed, the inner part of the accretion flow, where most of 
the luminosity is produced, shows the same density and temperature structure of the stationary models 
of spherical accretion onto black holes computed by Blondin (1986), Park (1990) and Nobili, Turolla & 
Zampieri (1991). Therefore, we can extrapolate the late-time evolution of the luminosity using their results. 
As shown in Figure [l^, the path described along the luminosity-accretion rate plane is in fair agreement 
with the solid line, which represents the analytic expression of I = l[fn) for stationary models derived by 
Blondin (1986) 

;~3xHr 7 f^) m 5 / 6 m>10. (53) 

Equation (^3|) is obtained assuming that all of the compressional work done by the gravitational field 
is converted into radiation in the inner accreting region where T oc r~ 5 / 8 (Blondin 1986). For rh > 10, 
Blondin's approximation provides a satisfactory fit to the numerical value of I = l(fn) calculated by Nobili, 
Turolla & Zampieri (1991). In Figure p~4| we plot the light curves of model IVa and IVb computed numerically 
along with their asymptotic behavior. For model IVa, the late-time accretion rate is estimated using 
equation (|50|). For model IVb, the late time rh is approximately calculated from rh(t) — rh(t re f)(t/t re f)~ 5 / 3 , 
where t re f — 6.5 x 10 3 i Qi o and rh(t re f) — 8.1 x 10 3 are taken from the computed model. Figure |lj shows 
that the analytic extrapolation is in fair agreement with the points computed numerically. 

In the interval 0.03 < rh <^ 10, Blondin's approximation is only approximately correct and is invalid 
below ~ 0.03, as the flow becomes transparent to its own radiation. To estimate the very late-time 
accretion luminosity one can adopt (as an order of magnitude estimate) the interpolation of the optically 
thin, stationary spherical accretion models (Park 1990; Nobili, Turolla, & Zampieri 1991) for which 

I oc to 2 rh < 0.03 . (54) 

In light of these findings, we can derive a simple scaling relation for the luminosity emitted by the black 
hole soon after the recombination front has propagated down to the inner accreting region. The scaling 
involves the main physical parameters of the flow at the onset of evolution. If the fluid is collisionless 
(hypothesis appropriate for model IVa) we can estimate the accretion rate according to equation (f30|) 



4tt 2 / 3 , ft 



-5/3 

m(t rec ) ~ p tock es I J . (55) 

Adopting equation (^3|) as scaling relation between / and rh, we find a luminosity at t rec 

-25/18 



1.2 x 10 



4 °(po<o) 5/6 ^J ergs" 1 . (56) 



Using as reference values those of Table |l| and a value of T rec ~ 10 4 K, for model IVa we estimate a 
luminosity L rec ~ 1.4 x 10 35 erg s _1 which is in close agreement we the numerical results. Equation (]5(j) 
thus provides an approximate expression for the accretion luminosity soon after clearing by recombination. 
L rec is a function of po, Tq and to at the onset of the explosion. Subsequently, the luminosity will decline 
with time as £ -25 / 18 (from equations pQ] and fl53]). 
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If the flow is radiation pressure dominated (model IVb) we can instead estimate the accretion rate at 
t ~ t rec from equation (BTh . Thus, we have 



ft \ _3 / 2 

m(t rec ) ~ Xi/2,GM bh ck es ^- I -Zp- ) (57) 



and the luminosity turns out to be 







(«)" 







-5/4 

L rec ~ 5.4 xl0 23 i £ 1^-1 (^J . (58) 

For model IVb L rec ~ 10 35 erg s _1 . The above relation can be applied to our numerical models since 
the density evolution is seen to preserve a self-similar character, even during the propagation of the 
recombination front: inside the accretion region, p(r,t) ~ po(t/to)~ 3 [r/r a (t)]~ 3 / 2 , where r a (t) ~ r a ,o(t/to) 
denotes the current value of the accretion radius. Equation ( |58| ) fixes the "level" of the luminosity at the 
moment of clearing of the envelope. Its time evolution will then depend on the nature of the flow after t rec . 
At late times we expect the transition to dust giving a luminosity L cx t~ 25 / 18 (as long as m ^ 10). 

In the next section we will use this scaling to give a first estimate of the time at which the bolometric 
light curve of SN1987A becomes dominated by the accreting black hole. 



8. Discussion and conclusions 

The computed light curves of model IVa and IVb show that, prior to recombination, the release of 
internal energy from the expanding stellar envelope dominates over accretion. Although during the first 
few hours the diffusive accretion luminosity can be very close to the Eddington limit, the diffusive and 
advection luminosities of the hot outer layers are largely superEddington. As shown by model IVb, during 
recombination the computed light curve shows the typical behavior observed in Type II— lb supernovae with 
low envelope mass. The luminosity at peak is comparable to that of SN1987A because the initial sound 
velocity of the hydrogen envelope is similar (see Tables and |). The high expansion velocity and low 
envelope mass of model IVb reduce all of the relevant timescales (see Table 0) so that recombination takes 
place at about 10 days after the onset of the explosion, instead of the 40-50 days inferred for SN1987A (see 
e.g. Woosley 1988). Nevertheless, our numerical results clearly indicate that no distinct signature in the 
light curve is found that could reveal the accreting black hole during this initial evolutionary phase. Only 
after the sharp decrease in luminosity due the clearing of the envelope by recombination, the accreting 
black hole starts to power the bolometric light curve of the computed models and I decreases as t~ 25 / 18 (see 
equations JSQ] and [j53| ). 

Slightly after maximum, SN1987A was powered by the radioactive decay of 56 Co and 57 Co and, at 
present, its light curve is consistent with emission from the decay of 44 Ti (Suntzeff 1997; see Figure |l4|). Can 
we discern the luminosity emitted by the accreting black hole above the contribution resulting from 44 Ti ? 
As already noted, the accretion history is very sensitive to the value of the sound speed (see equation |52| ) . 
Depending on the actual value of c Sj q in the inner part of the ejecta surrounding the compact remnant (the 
mantle), the accretion rate and, in particular, the total accreted mass vary significantly. An upper limit to 
the accretion luminosity can be inferred taking the initial parameters of the mantle quoted by Chevalier 
(1989; see Table |^). In this case, c Si o — 3 x 10 7 cm s _1 and t a> o ~ to, giving substantial fallback. The mass 
accretion rate can be computed using equation ( |50"| ) and the luminosity according to equation (|53|). Thus, 
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we have 

, . x -25/18 

l~3xlO- 7 ( Po t ck es ) 5 / 6 (-) . (59) 

Following Chevalier (1989), we adopt po^o ~ 10 9 C S S an d to ~ 7000 seconds, yielding a maximum accretion 
luminosity 

8 x 10 -3 

1 ~ (A/ Wl /Af Q ) 1 /3^( y ears)] 25 / 18 " (60) 

Equation ( |60| ) implies L~ 5 x 10 34 ergs~ 1 after t « 10 years, which is well below the present day bolometric 
luminosity of the remnant (~ 10 36 ergs _1 ; Suntzeff 1997) and also smaller than the luminosity estimated to 
result from radioactive decay. Thus, there is no observation that rules out the possibility that a black hole 
resides inside the SN1987A remnant. At very late times, when m <J 0.03, the corresponding upper bound 
would be 

/ t \ ~ 10/3 

l~2x 10- 8 ( Po t ck es ) 2 f-J , (61) 

or, using Table [|, 

7 x 10" 8 

~ [t(10 3 years)] 10 /3- ( > 

In Figure [l4] we plot the late-time light curve of SN1987A calculated using equations (60) and (p2|). As 
radioactive decay plummets at around 270-2700 years, it is clear that the black hole would only appear 
after about 900 years irrespective of the detailed numbers. After this time has elapsed, the luminosity of 
the remnant would be ~ 10 32 ergs _1 , too dim to be detectable using present technology, and possibly even 
a challenge for our distant descendants, as it would be hard to distinguish from the multitude of low mass 
stars crowding its field in the Large Magellanic Cloud. 

Since we have carried out a frequency-integrated calculation, we can only estimate the emission 
properties from the effective temperature and the photospheric gas temperature. During the optically thick 
accretion phase (when m > 0.03), the numerical results show that T e // and T p h are very close and have 
a value ~ 8000-10000 K. This corresponds to a mean photon energy in the few cV range, for which the 
bulk of the emission is expected in the visible band. Since T e ff and T p h are close together the continuum 
spectrum should be roughly a blackbody. This is not surprising because the emission processes are thermal 
and the gas in the inner accreting region is optically thick. Only at very late times [t ^> 2500 yr), the 
accretion flow becomes optically thin and the gas temperature reaches very high values (10 8 -10 9 if ; this 
regime has not been treated in the present numerical computation). So, at the time of emergence, the inner 
part of the accretion flow will be still optically thick and will be emitting roughly a black body spectrum 
peaked in the visible band. 

A number of interesting issues can be addressed within the framework presented in this paper. Can the 
presence of a stellar black hole in SN1987A be deduced from the excess accretion luminosity it produces ? 
What is the observational signature of the explosion of more massive progenitors (those with M ~ 30M©) 
on the light curve of the remnant ? Woosley, & Timmes (1996) suggested that a distinguishing signature 
would be a bright plateau that steeply falls to very low or zero luminosity, similarly to what found in our 
models. 

A number of assumptions limit the applicability of our results, amongst, the hypothesis of spherical 
symmetry and the simple chemical composition considered. As a result of an asymmetry in the explosion 
mechanism, black holes may come to birth with significant intrinsic velocities (a few x 100 km ~ 1 ). This 
is certainly observed for pulsars whose mean birth speed is around 250-300 km s _1 (Hansen & Phinncy 
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1997; Cordes & Chernoff 1997). However, since the initial sound speed of the ejecta can be of the order 
of 10 3 km s _1 , the accretion will proceed essentially subsonically. The central compact remnant tends 
to comove with gas expanding at its velocity and, if the ambient density is uniform, the conditions are 
independent of position within the core, as in the cosmological Hubble flow. Thus, even a relatively large 
intrinsic velocity should not strongly modify the results presented here. Deviations from spherical symmetry 
can arise if infalling matter has even a small amount of specific angular momentum, which presents a 
barrier to spherical infall (Chevalier 1996). In the case of SN1987A, the fact that the observed neutrino 
burst is consistent with implosion models without rotation (Burrows 1988) seems to indicate that angular 
momentum of matter initially close to the black hole was not very important (see also Chevalier 1996). In 
this case, the accretion should have proceeded almost radially, although mixing of outer mantle material 
with significant angular momentum into the inner region could modify our picture. 

A more realistic model with shell-like chemical composition and heating due to radioactive decay will 
enable us to address more quantitatively the problem of the visibility of black holes in otherwise successful 
supernova. This analysis will also set the basis for studying fallback onto a nascent neutron star in order to 
verify whether major outflows establish that reverse the process, during the decline of the supernova light 
curve. 
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A. Finite difference equations 

The full set of finite difference equations of relativistic hydrodynamics plus the radiation moment 
equations are summarized in this appendix. The indexes j and n will be used to denote spatial and time 
dependence, respectively. As discussed in the main text, p, B, wq and a are evaluated at mid-zones 
(j + 1/2), while r, M, u and w\ are evaluated at zone boundary (j). Mid-zone variables can be easily 
computed from zone boundary variables as: A/-1/2 = {Aj + A,_i)/2, where A is any variable. To compute 
zone boundary quantities from mid-zone quantities, it is necessary to take into account for the nonuniform 
spacing of the grid. Using linear interpolation, it is 

Aj = (4-1/2 + A j+1/2 )/(l + a) , (Al) 

where a is the fractional increment in grid spacing between successive zones (see Section ^). As far as time 
centering is concerned, u and w\ are evaluated at an intermediate time level (n + 1/2; time-shifted), while 
all the other variables are centered at the full time level (n). Linear interpolation/extrapolation in time has 
been used when necessary. 

Initial values for all the variables are specified on the grid according to what discussed in Section ||. 
Before starting the evolution, the radial component of the flow velocity u and the radiative flux w\ are 
advanced at time level i 1 / 2 solving equations (|^) and ( p4|) with a forward time integration scheme. 

The code starts to evolve first equation (||). Solving for the Eulerian radius at the new time level, we 
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obtain (j m in 5: i 5: jmax) 



„n+l 



At n+l/2 n+1/2 n+1/2 



(A2) 



where Ai™ +1 / 2 is the timestep between level n and n + 1 and 



n+1/2 



is interpolated at the zone boundary in 



space and extrapolated at time level n 
Equation 

(see e.g. May & White 1967; j mm + 1 < j < j 



1/2. Then, the time evolution of the matter density is computed. 



can be formally integrated in time and solved for P™^l/ 2 usm S the Crank-Nicholson operator 



Pi-l/2 - 



— n 
1/2 - Pj-1/2 



(r 2 )]_ 1/2 l - f p /2 

-F P ir 



where 



At 



1+1/2 n+1/2 
a 3 -l/2 



n+1/2 n+1/2 



3 



"j-1 



1+1/2 



-1/2 



2tt^ 



+ 1/2 



Ml 



-1/2 



n+1/2 

r i-i 



{wi) r - 



+ 1/2 

-1 



r 



n+1/2 
3-1/2 



(A3) 



(A4) 



In equation (A3), (r 2 )™^^ nas been computed in order to give the correct total volume of shell j — 1/2 



(May & White 1967). To derive equation (A3) and (A4), we used equation (|5|) and the two following 
relations, valid at constant Lagrangian time t: 



d_ 

dfi, 



bT 

d_ 



(A5) 
(A6) 



In equation 



n+1/2 , pu- 
tt. , ,„ and 1 . 

3-1/2 3- 



-1/2 
1/2 



have been extrapolated at time level n + 1/2. 



After calculating the optical depth r™ +1 and the Eddington factor with the advanced value of the 
matter density, the energy and the 0-th moment equations (equations and ^6|) are solved for B*^, 2 
and {w Q )™+l /2 (j 

rain ~t~ 1 ^ 3 — Jmax ) 



n+1 



-1/2 c J-l/2 
At n+l/2 
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(A7) 
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whcre B^l'il = {B?_ 1/2 + B]+l /2 )/2 and [w^l'/l = [(u*>)?_ 1/2 + {w )]^, 2 }/2. Since e, p and k P are in 



general rather complicated functions of p and B (see e.g. equations KM, |L5|]), equations ( |A7| ) and (A.8) 
form a highly non-linear system that has been solved by means of the Newton-Raphson method. To write 
the 0-th moment equation in finite difference form, r t t has been substituted with au (equation S) and the 
term b t t/b = — (pr 2 ).t/pr 2 has been expressed using the continuity equation (|j). 

The time slice at constant t is a characteristic direction for equation (|^) so that it can be formally 
integrated along the grid at the new time level t n+1 . Using the Leith-Hardy operator to approximate the 
exponential exp(— F a ) — exp{— /[(e,^ — hp^ — bs\) /hp]dp} in the interval [pj-1/2, Pj+1/2] (see May & 
White 1967), it reads (j min + 1 < j < j max ) 



l j-l/2 



h n +l /2 l + F a + F 2 /2' 



(A9) 



where 



F„ 



(hp) 3 



P3 ( e j+l/2 - £3-1/2) - — (Pj+1/2 - Pj-l/2) 
Pj 



Pj+1/2 ~ Pj-l/2 



s i)j 



n+l 



(A10) 



In equation ( A10| ), pj, pj and (hp)j are linearly interpolated at the zone boundaries. The difference 
Pj+1/2 — Pj-1/2 can be written as (1 + a)Apj_ 1 / 2 /2. Finally, the value of wi in s± has been linearly 
extrapolated at the new time level t n+1 . The boundary condition at p = Pj maw fixes a Jmax = 1 (equation 



After that the new timestep Ai™ +3 / 2 has been computed, the radial component of the fluid 4-velocity 
is calculated from the Euler equation (S) (j m in + 1 < 3 : < jmax — 1) 



™+3/2 = n+l/2 At n+l a n+l J ^ 



47rr„- 



+ ( 3 +fj ) Oo)j 



87rr j ( P3+1/2 - P3-1/2 
1 + a V Apj-1/2 

\ n+l 



(All) 



where Af n+1 is defined below (equation |A19|) and Tj, fty, pj, pj and /j have been radially interpolated 



using equation (Al). Along with wi, also the effective gravitational mass M has been linearly extrapolated 



in time at level t n+1 . At this point u" +1 can be interpolated from u" +1 ^ 2 and u™ +3 ^ 2 and the value of 
can be updated. The inner and outer boundary conditions for u are (equations p9[ and |p0||) 



1+3/2 _ n+l/2 . „+! „+! lu j min 
"jmin - U Jmin + ^ l a jmin f^n+1 ^2 



n+3/2 n+l/2 

u,- =u- 

Jmax Jmax 



(A12) 
(A13) 



where a" +1 has been radially extrapolated from the neighboring mesh points. 



The first moment of the specific intensity is computed from equation (|24|). In finite difference form it 
reads (j mm + 1 <j < jmax - 1) 
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where 
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(A16) 



As in the 0-th moment equation, to derive equation (A14) we have expressed r t and b^/b = — {pr 2 \ t /pT 2 
using equations (^|) and (^). To evaluate the radiation self-gravity term (the first term in curly brackets in 
equation [ A16fl ) an extrapolated value (in time) of wi is used. The inner boundary condition is a 5 point 
Lagrangian extrapolation in radius. Following Mihalas & Mihalas (1984) and Shapiro (1996), we impose 
the outer boundary condition applying equation (24) across the half-zone from j m ax ~ 1/2 to j max and 
substituting (wi)j max — g(wo)j max for (wo)j max in the gradients with respect to /i. The closure boundary 
factor g is given by equation (|l3|). 

Finally, the last quantity to be computed is the gravitational mass M. Integrating equation (^) along 
the time slice at constant t, we obtain (j m in + 1 < j < jmax) 



M™ +1 = M^l 



r n+l _ n+l 



Pj-1/2 (l + e :/-l/2) + (W ) j-1/2 + 
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The boundary condition at /1 = fij min fixes Mj min = Mq (equation p8| ]; see main text for explanation). 
The timestep is computed at each cycle before the integration of the Euler equation according to 



Ar+3 /2 = min 



At p ,At B ,At c ,(1.2At n+1 / 2 ) 
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where (j rnin + 1 < j < j max ) 



At p = min 



Ats = min 



0.05 



0.1 



n+l 

Pj-1/2 



n+l 

Pj-1/2 



_ M n+l/2 



P'- 



3-1/2 



nn+1 
3-1/2 



nn+1 
B j-l/2 



B j-l/2 



_ M n+l/2 



At r 



0.4- 



A A*j-i/2 



Courant condition . 
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Here v c = (f + 1/3) 1 / 2 is the radiation characteristic speed (the fastest speed on the grid). 

Artificial viscosity is inserted into the finite difference equations adding a scalar stress q j-1/2 to the gas 
pressure Pj-1/2 (von Neumann & Richtmyer 1950): 
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Fig. 1. Light curve for model I. The observed luminosity / = L/LEdd is plotted versus time t in units of 
the initial accretion timescale t a fi. The squares are the points computed numerically, while the solid line 
represents the analytic solution by Arnctt (1980). The arrow marks the diffusion timescale tdifj- 

Fig. 2. — Light curve for model II. The conventions and symbols are the same as in Figure 1. 

Fig. 3. — Evolution of the accretion luminosity and accretion rate at a fixed radius ro = 0.3r aj o, where 
r a ,o = GMbh/cl is the initial accretion radius, for model II. (a) Luminosity i(ro) versus time t. (b) 
Accretion rate rh = Mo/MEdd versus time t. 

Fig. 4. — Properties of the flow at selected times for model II. (a) Radial component of the 4-velocity |u| 
in units of the initial sound velocity c s ,o versus radius r in units of r a ,o- The dashed line denotes negative 
values, (b) Density p in units of the initial density po versus radius, (c) Gas temperature T versus radius, 
(d) Luminosity / versus r/r 0j o- At the onset of evolution the radiative flux is taken to be zero. Curves are 
labeled by time in unit of t a .a- 

Fig. 5. — Light curve for model III. The conventions and symbols are the same as in Figure 1. 
Fig. 6. — Luminosity / versus radius r/r a . n at selected times for model III. 

Fig. 7. — Computed light curve for model IVa. The observed luminosity I = L/L E dd (squares) is plotted 
versus time t in units of the initial accretion timescale t a a . When recombination takes place (i/i o ~ 42) 
the internal energy is released and the luminosity suddenly increases. After recombination the luminosity 
falls off exponentially until the emission is dominated by accretion (note the slow power law decline in I al 
late times). 

Fig. 8. — Properties of the flow at selected times for model IVa. (a) Gas temperature T versus radius. 
Note the formation of the recombination front (t/i a ,o ~ 42) and its rapid inward motion. The front stalls 
as soon as it approaches the innermost accreting region, (b) Luminosity I versus r/r a fi. (c) Radius of the 
recombination front r rec as a function of t/t a ,o- 

Fig. 9. Light curve, effective temperature and photospheric temperature for model IVb. (a) Luminosity 
I (squares) versus time, (b) Effective temperature T e ff (squares), as defined by equation (48), and gas 
temperature at the photospheric radius T p h (circles), as defined by equation (49), versus time. 

Fig. 10. — Early evolution of model IVb. (a) Radial component of the 4-velocity |u| versus radius r. The 
dashed line denotes negative values. Note the sign inversion of the velocity of the intermediate shells caused 
by the radiative force, (b) Luminosity I versus radius. The strong transient superEddington flux produced 
early in the evolution (t ~ t a> o) pushes outwards the marginally bound shells near the accretion radius and 
then propagates outwards. 

Fig. 11. — Accretion rate m(ro) at a fixed radius (r = 0.1r 0i o) versus time for model IVa (squares) and 
model IVb (circles). 

Fig. 12. — Mass flux rh versus radius at t = 6 x 10 3 i a ,o for model IVb. Note the constancy of rh with radius 
in the inner accreting region. 

Fig. 13. — Tracks on the luminosity-accretion rate plane (l-m) for model IVa and IVb. The observed 
luminosity I is plotted as a function of the accretion rate rh (evaluated in the inner accreting region). 
Squares denote the computed points for model IVa, while circles denote those of model IVb. The solid line 
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is the analytic expression of I = l{m) for stationary models derived by Blondin (1986; see equation (52) in 
the main text). 

Fig. 14. — Computed and observed bolometric light curves for SN1987A. The luminosity L is plotted as 
a function of time t. Squares {circles) give the solution for model IVa (IVb). The dotted lines represent 
the extrapolation of the late-time evolution for the computed models (see Section 7). The triangles are the 
bolometric luminosity of SN1987A (data taken by McCray 1993, Woosley & Timmes 1996, Arnett 1996 and 
Suntzcff 1997). The dashed lines represent the expected contribution from the decay of radioactive elements 
(0.07 M of 56 Co and ~ 5x 10~ 5 Af Q of 44 Ti). Finally, the upper dotted line denotes the expected bolometric 
luminosity emitted by a putative black hole in SN1987A. The arrow marks the time (t ~ 900 years) at which 
the upper dotted line crosses the dashed line (radioactive emission from 44 Ti): from this moment on, the 
bolometric luminosity of SN1987A would be dominated by the energy released by accretion. 
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TABLE 1 

Parameters for the computed models 



model fc a c s , (cms- 1 ) M cloud (M ) r™, (cm) Vo b ( cms_1 ) p c (g cm" 3 ) T d (K) t a>0 e t f t d ,// s ttran S h tree 1 



I 1 3.5 X 10 7 10~ 5 3 X 10 12 6.4 x 10 s 1.8 X 10~ 10 9 X 10 4 1.3 h 1.3 h 1.6 h 19 h oo 

II 1 3.5 X 10 7 10~ 2 3 X 10 12 6.4 x 10 s 1.8 X 10~ 7 5 X 10 5 1.3 h 1.3 h 2.1 d 25 d oo 

III 0.1 3.5 X 10 7 10~ 2 3 X 10 12 6.4 X 10 7 1.8 X 10~ 7 5 X 10 5 1.3 h 13 h 6.7 d 250 d oo 
IVa 1 3.5 X 10 7 10~ 2 3 X 10 12 6.4 x 10 s 1.8 X 10~ 7 5 X 10 5 1.3 h 1.3 h 2.1 d 25 d 2.7 d 
IVb 0.1 10 s 1 3 X 10 12 1.5X10 9 1.8 XlO- 5 2.7 X 10 6 0.055 h 0.55 h 14 d 107 d 6.2 d 



b Expansion velocity of the outermost shell. 

c Density of the cloud at the onset of homologous expansion. 

d Initial temperature in the inner part of the cloud. 

e Accretion timescale (equation [33]). 

f Expansion timescale (equation [34]). 

s Diffusion timescale (equation [35]). 

h Thick— thin transition timescale (equation [36]). 

Recombination timescale (equation [37]). 

Note. — The mass of the black hole is taken to be = 1.5Mq. 



1 



TABLE 2 



Parameters of SN1987A 



fc a c s , (cms- 1 ) M cloud (M Q ) r out (cm) Vo b ( cms_1 ) Po c (g cm" 3 ) T d (K) t a>0 e t f t d ,// S *tran S h tree' 



Parameters for the hydrogen envelope (taken from Arnett 1996) 



0.0067 


1.6 X 10 8 


15 


3 X 10 12 


4 X 10 8 2.7 X 10~ 4 6.7 X 10 6 


0.014 h 


2.1 h 


103 d 


4.2 y 


58 d 








Parameters 


for the mantle (taken from Chevalier 1989) 












0.66 


3.5 X 10 7 


4 


9 X 10 11 


1 X 10 8 3 X 10~ 3 5.7 X 10 6 


1.3 h 


2 h 


107 d 


6.2 y 


110 d 
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